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ABSTRACT 

We present a general formulation for stability analyses of radiative shocks with multi- 
ple cooling processes, longitudinal and transverse perturbations, and unequal electron 
and ion temperatures. Using the accretion shocks of magnetic cataclysmic variables 
as an illustrative application, we investigate the shock instabilities by examining the 
eigenfunctions of the perturbed hydrodynamic variables. We also investigate the effects 
of varying the condition at the lower boundary of the post-shock flow from a zero- 
velocity fixed wall to several alternative type of boundaries involving the perturbed 
hydrodynamic variables, and the variations of the emission from the post-shock flow 
under different modes of oscillations. We found that the stability properties for flow 
with a stationary-wall lower boundary are not significantly affected by perturbing the 
lower boundary condition, and they are determined mainly by the energy-transport 
processes. Moreover, there is no obvious correlation between the amplitude or phase 
of the luminosity response and the stability properties of the system. Stability of the 
system can, however, be modified in the presence of transverse perturbation. The 
luminosity responses are also altered by transverse perturbation. 

Key words: accretion - - shock waves - - stars: binaries: close - - stars: white 
dwarfs 



1 INTRODUCTION 

The time-dependent properties of radiative shocks have been investigated by many researchers in different settings, for 
example the interactions between supernova shocks and the interstellar medium, and accretion flow onto compact objects 
(e.g. Falle 1975, 1981; Langer, Chanmugam & Shaviv 1981, 1982; Chevalier & Imamura 1982; Imamura, Wolff & Durisen 1984; 
Chanmugam, Langer & Shaviv 1985; Imamura 1985; Bertschinger 1986; Innes, Giddings & Falle 1987a, b; Gaetz, Edgar & 
Chevalier 1988; Wolff, Gardner & Wood 1989; Imamura & Wolff 1990; Houck & Chevalier 1992; Wu, Chanmugam & Shaviv 
1992; Toth & Draine 1993; Dgani & Soker 1994; Strickland & Blodin 1995; Wu et al. 1996; Imamura et al. 1996; Saxton, 
Wu & Pongracic 1997; Saxton et al. 1998; Hujeirat & Papaloizou 1998; Saxt on fc Wu 1999). Many of these shocks are found 



to be thermally unstable. For instance, a numerical study by Langer et al. (1981) showed that the post-shock accretion flow 
in magnetic cataclysmic variables (mCVs), binaries in which a magnetic white dwarf accretes material from a red-dwarf 
companion star, suffers thermal instabilities and hence fails to attain a steady-state. The accretion shock is driven to oscillate, 
giving rise to quasi-periodic oscillations in the optical luminosity. A similar conclusion was obtained by Chevalier & Imamura 



(1982), using a linear perturbative analysis. 

The stability of the radiative shock depends on the energy transport processes. In the case of accretion shocks in mCVs, 
thermal bremsstrahlung and cyclotron radiation are the most important cooling processes (e.g. Lamb & Masters 1979; King 
& Lasota 1979). Bremsstrahlung and cyclotron cooling have very different temperature dependences, and hence influence the 
stability properties of the accretion shock differently. 



In the stability analysis of Chevalier & Imamura ( 1982 ) , the total cooling effects were approximated by a single radiative 
loss term A oc p a T h depending on temperature T and density p. Various choices of the power-law indices (e.g. a = 0.5 and 
b = 2 for bremsstrahlung cooling) were investigated, and they have found that radiative shocks with larger b (i.e. stronger 
temperature dependence) are more stable against perturbations. 
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The individual cooling processes were s ubsequently considered explicitly in the numerical study of accretion onto magnetic 
white dwarfs by Chanmugam et al. (1985). An effective cooling term was constructed to mimic the effects due to optically 
thick cyclotron cooling. Their study showed that efficient cyclotron cooling stabilises the shock. Wu et al. (1992, 1996) further 
investigated the same system and found that in spite of the suppression of shock oscillations in the presence of cyclotron 
cooling, the oscillation frequency appears to increase quadratically with the magnetic-field strength. Moreover, each of the 
cooling processes, bremsstrahlung and cyclotron, dominates in about half of the phases of an oscillatory cycle, allowing the 
perpetuation of small-amplitude oscillations provided that the magnetic field is moderate or weak (B ^10 MG). 

Li near p erturbative analyses of accretion shocks with bremsstrahlung and cyclotron cooling were carried out by Saxton 
et al. ( 1997[ ). A composite cooling function (following Wu, Chanmugam & Shaviv 1994) was considered, comprising the 
sum of a term for b remsstrahlung cooling and an effective term for cyclotron cooling. The analysis was further extended by 
Saxton et al. ( 1998 ), in which the cooling function is a sum of terms for bremsstrahlung cooling (Ab r oc p 2 T 5 ) and a second 
power- law process with a destabilizing influence (A 2 oc p a T b for b > 1). The cases considered included that of a cooling term 
(A cy oc p 015r p 2 - 5 ^ effectively approximating the energy loss due to cyclotron radiation in the geometry and flow conditions of 
the post-shock regions of mCVs. It was found that a simple comparison of cooling and oscillation timescales is insufficient to 
understand the instabilities of the shock under various modes. 

When the radiative cooling is fast c ompared to the electron-ion energy exchange, the electron and ion temperatures are 
generally unequal. Imamura et al. (1996) considered bremsstrahlung and Compton cooling, and a more general perturbation 
of the shock in both longitudinal and transverse directions. Their study showed that the electron-ion exchange process and 
the presence of transverse perturbations can destabilize each mode compared to the purely longitu dinal and one-temperature 
cases. Saxton & Wu (1999) generalised the works of Chevalier & Imamura (1982), Imamura et al. (1996), Saxton et al. (1997) 
and Saxton et al. (1998) for radiative accretion shocks by considering multiple cooling processes explicitly, the two-temperature 
effects and transverse perturbations (as for a corrugated shock). In the case of mCVs the introduction of two-temperature 
effects complicated and broke down the strictly monotonic stabilisation of modes with increasing cyclotron efficiency in one- 
temperature shocks with bremsstrahlung and cyclotron cooling, and the influence of transverse perturbation was not always 
able to destabilize oscillatory modes in the presence of both bremsstrahlung and cyclotron cooling. (For a review of stability 
of accretion shocks, see Wu 2000.) 

The present paper expands upon the studies of the eigenvalue in Saxton & Wu (1999) to examine the amplitudes and 
phases of the eigenfunctions which describe the response of the post-shock structure to perturbations of the shock position. 
As an illustrative case we consider the accretion shocks of magnetic white dwarfs. Knowing the response of the hydrodynamic- 
variable profiles in turn provides information about other characteristics of thermally unstable shocks, such as the responses 
of post-shock emissions due to the shock oscillations. 



2 FORMULATION 
2.1 Hydrodynamics 

In accretion onto white dwarfs, the supersonic flow meets a stand-off shock where the inwards falling matter is abruptly 
decelerated to attain a subsonic speed. The shock sits above the white-dwarf surface at a height x s ~ ;j«ffi C ooi, where the 
free-fall velocity is vg = (2G , M w d/-R w d) 1/ ' 2 , the cooling timescale is t coo \ ~ n c k B T s /A, and A is a radiative cooling function. 
(M w( j and i? w d are the white-dwarf mass and radius respectively; fc B is the Boltzman constant; n c is the electron number 
density; and T s is the shock temperature.) 

The time-dependent mass continuity, momentum and energy equations for the post-shock accretion flow are 

(^+v-v) )O + p(V-v)=0, (1) 
9 +v-v)v + VP = 0, (2) 



dt 

| +v .v)p- 7 ^(|+v.v)p— (7-1JA, (3) 

(* + v.v)p.-7^(|+v.v)p=( 7 -l)(r-A), (4) 

where the total radiative cooling function A and the electron-ion energy exchange term F are local functions of the density p, 
total pressure P, the electron partial pressure P e , and the flow velocity v. The explicit form of the exchange function is 



4v / 27re 4 n c 7i i In C 



\3/2 



(5) 



where ni >e are the ion and electron number density, and #i jC = fc B Ti, e /mi jC c 2 with 71, c being the corresponding temperatures. 
The constants are the ion and electron masses, e is the electron charge, c is the speed of light and InC is the Coulomb 
logarithm (as in e.g. Melrose 1986). An adiabatic index 7 = 5/3 for an ideal gas is assumed, and the equation of state 
P = pfc B T//im H is considered, where m H is the mass of hydrogen atom. 



Stability analyses of two-temperature radiative shocks 3 



The total cooling function is written in a form of the bremsstrahlung-cooling term and a multiplicative term expressing 
the ratio of the losses due to the second cooling process and bremsstrahlung cooling. The second process is characterised by 
its power-law indices of density and electron pressure (a = b~ 1/2 and (3 = 3/2 — a + b for a general cooling term A2 oc p a T b ), 
and by the parameter e s , which is the relative efficiency evaluated at the shock. (Larger e s implies a more efficient second 
process.) 

A = A br + A 2 = A br [1 + e s /(r , tt c )] , (6) 
and a function is defined to relate the primary and secondary cooling processes: 



/(ro,7T e ) = 




where P CjS and p s are the shock values of electron pressure and density. The dimensionless parameter a s = (P c /Pi) s is the 
ratio of electron and ion pressures at the shock. The bremsstrahlung cooling term Ab r = Cp 2 (P c / p) 1 ^ 2 , where the constant is 
C = (27rfc B /3m e ) 1 / 2 (2 5 7re 6 /3/im e c 3 )(/it/fe B mp) 1 ^ 2 (7 B , with m p the proton mass, h the Planck constant, /1 the mean molecular 
weight and g B « 1 the Gaunt factor (see Rybicki & Lightman 1979). For completely ionised hydrogen plasma p — 0.5 and 
the constant has a value C w 3.9 x 10 16 in c.g.s. units. 



2.2 Perturbation 

A first-order perturbation is considered for the shock position x s and velocity v s : 

v s = v s ie aT , (8) 

x s = a; s o + Zsie " , (9) 

where to is the frequency, and k is the transverse wavenumber of perturbation in the y (transverse) direction. The shock is at 
rest in the stationary solution, « s o = 0, and the perturbed motion of the shock has v s i = x s iuj. The dimensionless frequency 
and transverse wavenumber are 

k = Xsok , (10) 
8 = ^cu . (11) 

VS 

The eigenfrequencies are complex, S — 5 R +iS I: with 5, being the dimensionless frequencies of the oscillations, and S K the 
stability term. When 8 R is positive the perturbation grows; when 8 R is negative the perturbation is damped. 

The post-shock position coordinate is labelled by £ ee x/x b , which is £ = 1 at the shock and £ = at the white-dwarf 
surface. The size of the perturbation is parameterised by e ee v s i/vg = Sx s i/x s o- The other scales to be eliminated are a; s o 
the stationary-state shock height and p a the mass density of the pre-shock accretion flow. The hydrodynamic variables are 
expressed as 

p{i,V,t) = Pa-Co(C)(l+eA c (C)e lfe! ' + " t ) , (12) 

v{H,y,t) = -v B -To(0((l+£X^y ky+ ^),sX y (^ k y + - t ) , (13) 

P{£,V,t) = p,v 2 B -M0^ + ^A0e ky+ ^) , (14) 

P.(Z,V,t) = Pa^-7r (C)(l+eA c (e)e lfc! ' + ^) . (15) 

where ^0, to, 7To and 7r c are dimensionless density, velocity, total pressure and electron pressure in the stationary solution; 
and A^, A T , A y , A^ and A c are complex functions representing the response of the downstream structure to perturbation of 
the shock height. These five functions describe perturbations of the density, longitudinal velocity, transverse velocity the total 
pressure and the electron pressure respectively. 

Separating the time-independent terms from the hydrodynamic equations yields two algebraic equations and two differ- 
ential equations for the stationary case: 

Co = 1/ro , (16) 

TTO = 1 - T(, , (17) 



d£ _ 77T0 — To 

dr A 



1 - Y ) ("fltO - T ) - 77T C 



(18) 
(19) 



dn e _ 1 

dn to \ A, 

The electron-ion energy exchange and cooling processes are described by appropriate dimensionless forms: 

f = (7 - lXpa^/Zso)-^ = ( 7 - l^cVci 1 '^^ - , (20) 
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A = (7 - l)(p a 4/x a o) _1 A = (7 - i)V» c 



[1 + e s / (r , 7T C 



(21) 



where the constant ^> c is determined by normalisation of the integrated stationary solution and the parameter f/> c i is a ratio 
between the electron-ion energy exchange and radiative cooling timescales, as described in Imamura et al. (1996) and Saxton 
& Wu (1999). The physical values of tp c and tp c i are given in Saxton & Wu (1999) and Saxton (1999). 
The first-order perturbation is determined by the matrix differential equation: 



dr 







X T 


1 


Xy 




X-jt 




. Ac . 





1 



-1 



1/ro 
-l/r 


1/jro 



-77ro/r 1 

— (77T0 — ro) /to 

7 —7 

7—7 7/V0 — (77To — To)/ToWe 

where the F functions, which are comprised of terms that do not include derivatives of the A variables, are given by 





" Fi ' 




F 2 




F 3 




Fi 




Fs 



Fi(ra,TT e ,(,) 
F 2 (to,ivc,0 

F 3 (t ,1Vc,0 
F4(To,7Te,£) 

^5(713,^6,0 



-£(lnTo) — 8X( + iKToXy , 

-{8 — tq)A t + ^(lnro)' + TqX^ — tq(A^ — A T ) 

-(8 — tq)X v + in(l — ro)A,r + inr^/S , 



-7r 5(A^ — 7Ac) - A 



+r 



TV e 5(X, 

5 



7A C ) 



A 

TO A 



3 1 

-ff2(-TO,7T c )A C + -£i(t ,7T c )A c 

3 1 

2ff2(r ,7r c )A c + -gi(T ,n c )X c 



X T 

A T - 



A* + 



A c + 



.1 

TO 
1" 
TO. 



27T C A C 



A T A c H - — 

T 



7T{) — 27T C 

where the primed quantities are derivatives in terms of £. The functions gi(ro,n c 
2e s a/(To,7r c ) 



and 32 (to, tv c ) are defined as 



ffl(T0,7T c ) = 1 + 



ff2(T ,7T c ) = 1 



1 + £ s /(t , 7T c 



e s /3/(r ,7r c ) 



1 + £s/(to, 7T c ) 



(22) 



(23) 
(24) 
(25) 

(26) 



(27) 

(28) 
(29) 



The complex matrix differential equation ( p2[ ) can be decomposed into ten first-order real differential equations in terms of the 
functions of the stationary solution f (to) and 7r e (To), and the real and imaginary parts of each of the perturbed variables. It 
can be shown that equation (^) is the general description and that the more restricted formulations in Chevalier & Imamura 
(1982), Saxton et al. (1997, 1998) can be recovered from it under specific assumptions (Appendix |a|) . 



2.3 Cooling functions 

Optically-thick cyclotron cooling depends not only on the local properties of the hydrodynamic variables but also on the 
radiation field. A self-consistent description of cyclotron loss usually requires solving the equations for radiative transfer and 
the hydrodynamics simultaneously. The particular geometry and physical conditions of the magnetically channelled accretion 
flow in mCVs, however, permit a simplification (see e.g. Cropper et al. 1999). A functional fit involving the density and 
temperature dependences of the cut-off frequency yields an approximate cooling term A cy oc p 015 T e 2 ' 5 (see Langer et al. 
1982; Wu et al. 1994). 

In our analysis the effective cyclotron cooling term has power-law indices (a, 0) = (2.0, 3.85), and the efficiency parameter 
e s (Wu 1994) depends upon the temperature, density, magnetic field, and geometry of the emission region. In the limit of a 
one-temperature accretion flow the efficiency of electron-ion energy exchange ij) c i is large, and the ratio of pressures a s tends 
to unity. For two-temperature shocks, a s is determined by the electrons carrying energy into the region above the ion shock, 
which is a complication beyond the scope of this paper. Following Imamura et al. (1996) and Saxton & Wu (1999), we treat 
a E as a parameter. 



2.4 Stationary structure 

We assume the stationary-wall condition, which requires a zero terminal velocity at the lower boundary. The stationary 
solution can be obtained by direct integration, after substituting equations (16) and (17) into equations (18) and (19). In 
Figure |l|we show the stationary velocity structures (ro) of the post-shock region for various choices of the system parameters. 
The electron and ion sound speeds (c c , q) are also plotted on the same scale. The density is related to the velocity by 
po = Ph/to and the respective temperatures are proportional to the squared sound speeds. Two-temperature effects are more 
significant in cases when cyclotron cooling is efficient (e.g. e s — 100), and when the electron-ion exchange is inefficient (i.e. 
small ip e i). When two-temperature effects are unimportant, we recover the velocity, density and temperature structure of the 
one-temperature calculations (Wu 1994; Chevalier & Imamura 1982). A detailed discussion on the two-temperature stationary 
structures of the post-shock flows and their emission will be presented elsewhere (Saxton, Wu & Cropper, in preparation). 
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Figure 1. Statinary structures of the post-shock flow, with to the flow velocity normalised to vg (dashed lines) as functions of the 
normalised position £. The electron and ion sound speeds are represented by the upper and lower solid curves respectively. The upper, 
middle and lower rows represent parameter choices (<7 s ,V> c i) = (0-2,0.1), (<r s ,V , ei) = (0.5,0.5) and, (<r s ,i/ , oi) = (1> 10) respectively. The 
columns represent values e B = 0, 0.1, 1, 10, 100 from left to right. 



2.5 Time-dependent solutions 

For the time-dependent solution, the set of differential equations for the perturbed variables are integrated numerically from 
the shock down to the lower boundary for trial values of the complex eigenfrequency 8. The variables at the shock are 
determined by the shock-jump conditions: A^ = 0, \ T = —3, \ y = 3in/S, and = A e = 2. The condition at the lower 
boundary is less well defined. The stationary-wall condition only requires a zero terminal velocity at the bottom for the 
stationary solution. A lower boundary condition for the perturbed variables, such as the requirement that the flow should 
stagnate (A T = 0), the total pressure be constant (A w = 0), or some other condition relating density, pressure and longitudinal 
velocity (e.g. A ? + A T =0), may be chosen. 

The time-dependent solution can be expressed as a linear combination of oscillations of different eigenmodes. For a 
hydrodynamic variable X, we have 

X(£,t) = Xoi^Y, [ a « A *."(0exp (^<w) 

n 

where Xa(t£) is the stationary solution, a„ is the relative strength of the mode n, 8 n the eigenfrequency and A x , ra the eigenfunc- 
tion. While the complex-<5 eigenvalues provide information about the global stability properties and oscillating frequencies of 
the modes, the A-eigenfunctions describe the more local dynamical properties. The absolute value of the A-function determines 
the local oscillation amplitude of the mode, and the phase indicates whether the oscillation of the hydrodynamic variables 
lags or leads in one region with respect to another. 



(30) 
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3 EIGENVALUES 

The eigenvalues for two-temperature oscillating shocks with a "perfect" stationary-wall lower boundary condition, i.e. r = 
A T = (see Saxton 2001), have been discussed in our previous paper (Saxton & Wu 1999, see also Imamura et al. 1996). More 
restricted studies on the one-temperature radiative shocks were presented in Saxton et al. (1997) and Saxton & Wu (1998) 
(see also Chevalier & Imamura 1982). Particular results of these studies include 

(i) In the one-temperature case, the frequencies are quantised like modes of a pipe open at one end, 5 1 m 8 IO (n — 1/2) + 8 C 
with a small correction S c . When two-temperature effects are strong the "stationary- wall" condition at the lower boundary 
loses importance, and the frequency quantisation becomes more like that of a doubly-open pipe, i.e. <5 : m S IQ n. 

(ii) The frequency spacing S IO decreases as the efficiency of cyclotron cooling (e s ) increases, but tends to increase when the 
efficiency of electron-ion energy exchange (t[>ei) decreases. 

(iii) Increasing e s generally stabilises the modes, but when two-temperature effects are extreme there are situations in which 
an increase of e s destablises some modes. 

(iv) In presence of transverse perturbation, there are maxima of instability at certain values of transverse wavenumber k 
for each longitudinal mode. However, there are values of (a s ,ip e i, e s , a, /?) for which some longitudinal modes are stable at all 
k ranges. 

(v) For a given longitudinal mode, the shocks are generally stable against transverse perturbation of sufficiently large k. 

(vi) Two-temperature effects affect the stability of cyclotron-cooling dominated shocks (e s ^> 1), but have less influence 
when bremsstrahlung cooling dominates (e s ~ 0). 



4 EIGENFUNCTIONS 
4.1 Amplitude-profile 

In one-temperature accretion shocks, the non-trivial eigenfunctions are A T , A-„-, and A^, corresponding to total pressure, 
longitudinal velocity and density. In the two-temperature shocks, the non-trivial eigenfunctions also include that of the 
electron pressure A e When there is a transverse perturbation, a transverse perturbed velocity eigenfunction X y becomes 
important also (see §(?]). 

Our choice of boundary conditions requires that the perturbed longitudinal velocity A T has a fixed value of —3 at the shock, 
and a stagnant-flow at the white-dwarf surface implies a zero A T at £ = 0. (Here and hereafter in this section, unless specified, 
the stagnant-flow condition is assumed.) The amplitude profiles generally have local minima and maxima (see Appendix B), 
which are analogous to the nodes and antinodes in the oscillations of a pipe. The number of nodes and antinodes depend on 
the harmonic number n, and their positions are determined by the radiative-transport processes and the system geometry. 
Unlike the nodes of an ideal pipe, the amplitude minima of oscillating shocks do not always go to zero, and the amplitude 
maxima do not have the same values. Moreover, the nodal features of one hydrodynamic variable may not coincide with those 
of another hydrodynamic variable. 



4-1.1 \\ T \-profile 



The | A T | profiles do not show very distinctive nodes and antinodes. The amplitude |A T | increases very rapidly with £ from a 
zero value at the lower boundary, reaches a plateau or a peak at larger £ and then declines with a gradual tail towards the 
shock boundary value. The rapid rise of |A T | occurs only in a small region where £ ^ 10 -4 , and the size of the region seems to 
be independent of the harmonic number n (see Figure El) . The profiles for e s = and 1 are similar, with the plateau value of 
At | of the fundamental mode « 3 (e.g. Figure |B5| ). For higher harmonics, the plateau values of |A T j are larger. The maxima 
(antinodes) generally become more visible when n increases (e.g. Figures B6, B7, B8). 

For the cases with large e s , |A T | has a strong peak near the lower boundary. The peak heights tends to increase with n. It 
is worth noting that in the extreme of small tp e i, large <r s and large e s (strong two-temperature effects and equal temperature 
for electrons and ions at the shock), a small peak in |A T | is also present near the shock (see Appendix D.2.4 in Saxton 1999). 



4-. 1.2 \\Tr\-profile 

The boundary value of is 2 at the shock, but it is not defined at the white-dwarf surface for the stagnant-flow condition. 
The |A^| profiles do not show distinctive nodes and antinodes fo r n = 1 and 2. The nodes and antinodes, however, become 
more visible when n increases and ip e i becomes small (Figure 



4.1.3 A c I -profile 



The electron pressure eigenfunctions A e show nodes and antinodes. The sharpness of the nodes depends on the strength of the 
electron-ion energy exchange, t/'ei- When the electron-ion energy exchange is efficient, the nodes are weak. When the electron- 
ion energy exchange is inefficient, (e.g. _^ci = 0.1), the amplitude minima are narrow in £ and deep in terms of amplitude. 



(See e.g. the e B = 100 curves of Figures B1-B4) 
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Figure 2. The amplified view |A T | profiles (left column) and the phases ip T (right column) near the lower boundary for the case with 
(cr B ,ip e i,e B ) = (0.5,0.5,1.0) and a stagnant-flow boundary condition (A T = at £ = 0). The |A T | profiles and the phases of the modes 
n = 1, 2, 3 and 4 are shown in rows from top to bottom. The |A T | profiles are similar for the four modes, but the phases are very different. 



For two-temperature shocks, it is the electron pressure ir c rather than the total pressure no that appears in the cooling 
terms and electron-ion exchange term. The electron pressure eigenfunction A c therefore adopts the role that is played by 
the total pressure eigenfunction \ n in the one-temperature case. When the electron-ion energy exchange is weak (i.e. small 
ipei), the A w and A c eigenfunctions are very different. The different modes have sharp and strikingly different features in A c |, 
whereas the total pressure eigenfunction is relatively featureless and smoothly varying. When the energy exchange between 
electrons and ions is efficient (i.e. large ipei), the disparties between their pressures becomes small throughout most of the 
post-shock region. When the system approaches the one-temperature condition, and the total pressure follows the behaviour of 
the electron pressure, and the variations of the electron pressure and total pressure eigenfunctions, A c and \ n , are comparable 
in amplitude over all £. 

Because of the lower boundary conditions that we assume, the electron and ion temperatures are zero and the electron and 
ion pressures always become equal at £ = 0. The relative oscillation of the electron pressure must equal the relative oscillation 
of the total pressure at the lower boundary. I.e. A e (£) — > A,r(£) as £ — > 0. However, there is no condition to determine the 
particular value at which these quantities must meet at £ = 0, for any given mode under particular conditions. 



4.I.4 \X ( \-profile 

The density eigenfunctions, A^, have strong and distinct nodes and antinodes. For each mode the positions and magnitudes 
of the node-antinode-node sections are affected by the efficiency of cyclotron cooling e s . When e s is small (£, 1), the nodes are 
nearly evenly spaced and the antinodes near the shock tend to have lower amplitudes than those nearer to the white-dwarf 
surface. When e s is larger, the nodes near the shock are more widely spaced and the antinodes near the shock tend to be higher 
in amplitude. The innermost antinode in the region £ <, 0.1 becomes higher as e s increases through the cases of e B = 0, 1, 100 
studied here. 

Because of the sharp node and antinode features appearing in the profiles of |Af | (and |A e |) we suspect that the density 
and electron pressure may be the quantities which determine the oscillatory properties of the shock. These hydrodynamic 
variables appear explicitly in the functions for the cooling and electron-ion energy exchange processes. Other hydrodynamic 
variables (i.e. the flow velocities and the total pressure) have amplitudes varying less rapidly in £, probably indicating a less 
active involvement in determining the oscillatory behaviour of the shock. 

For the studied cases of (cr s , tp e i, £s), the modes divide into two classes based on the qualitative features of their amplitude 
profiles. Profiles of the fundamental and first overtone (n = 1, 2) are more alike than the profiles of higher modes. The n > 2 
modes have n — 2 nodes at intermediate £ positions, whereas the n — 1, 2 profiles are slowly-varying in £. The actual number 
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of density nodes in each profile is n — 1, including the shock (where A^ = as a boundary condition). An additional shallow, 
and usually indistinct, node occurs near the fixed wall boundary (small f) for some cases of large e s . When two-temperature 
effects are strong, |A e | also has sharp nodes for n > 2 modes. 



4-1.5 nodes, antinodes and sound speeds 

The dimensionless electron and ion sound speeds are given by 
and 



(31) 



Ci = ^/7 r o(l - T - 7T C ) (32) 

respectively, and the mean sound speed is 

Cs = \J"1T {1 - T ) . (33) 

The sound speeds determine the local dynamical time and length scales, which influence the mode spacing. In regions of the 
post-shock flow where the sound propagation is fast, nodes of the eigenfunctions tend to be further apart. 

The sound speeds generally decrease from the shock down to the lower boundary (see Figure |l]). When the gradients of the 
sound speed have less dramatic change (e.g. the case e s = as in leftmost column of Figure |l|) the nodes of the eigenfunctions 
are more evenly spaced. When e s is large the sound speeds are proportionately larger in the upper region of the flow, resulting 
in wider node spacing in the high-£ region (contrast the different e s curves of |A^| in Figure |B4| ). 

The local amplitudes of antinodes of the density eigenfunctions are greatest in regions near the shock for large e s . This is 
counterintuitive because previous studies (e.g. numerical simulations by Chanmugam et al. 1985, and Wu et al. 1992, 1996) 
show that oscillations are globally suppressed by cyclotron cooling. This global result would lead us to expect lower local 
amplitudes in the cyclotron-dominated region, and reduced amplitudes in systems with greater e s . We believe that there is a 
connection between the increased antinode amplitudes and the increased node spacing in the cyclotron-dominated region when 
e s is great. The amplitude enhancement is not related to two-temperature effects, because it occurs in the one-temperature 
extreme (large ip e i and a s — > I) as well as the general two-temperature cases. The behaviours of the nodes and antinodes (of 
the density eigenfunction) probably depend more sensitively on the overall form of the stationary solution than on the energy 
exchange processes present in each region of the flow. 



4.2 Phase 



The phases of the various perturbed hydrodynamic variables are given the labels ip^, ip T , ip y , ip n and <p e for density, longitudinal 
velocity, transverse velocity, total pressure and electron pressure respectively. The zero values correspond to oscillations in 
phase with that of the shock height. 



4-2.1 Boundary conditions and phase 

At the shock (£ = I) there are boundary conditions on all of the A variables. The only phase which is not explicitly determined 
is tpQ. Our calculations, however, show that the density oscillation is approximately in quarter-phase ahead of or behind the 



shock height oscillation, except when e a ~ 100 (see Figures B2, B6) 



There are fewer boundary conditions on the hydrodynamic variables at the bottom of the post-shock region. For the 
stationary- wall condition, only the longitudinal velocity variable is constrained, |A T | — > as £ — > 0. The phase tp T is not 
determined. The other perturbed variables have no definite lower boundary values, and their phases depend completely on 
the energy transport processes. 

The phases ip c and tp„ are not completely independent. Since the electrons and ions have reached the same temperature 
at the lower boundary, the oscillations of electron pressure and the oscillations of total pressure are identical at £ = 0. The 
variables X n and A c have the same values at the shock and the base but different values in between. 



4-2.2 Interpretation of dip / d£ 

If the phase increases with £, we consider the oscillation to be propagating downwards from the shock towards the lower 
boundary. (This is because the oscillations in the high-£ region lead the oscillation in the lower-£ regions.) If the phase 
decreases with £ then the oscillation is considered to be propagating upwards. The upward and downward propagation are 
labelled positive [+] and negative [— ] respectively. If a complex A-function is viewed along the £ axis, with the real and 
imaginary A parts horizontal and vertical, then the complex function appears to wind about the origin in either a clockwise 
manner for positive propagation, or an anticlockwise manner for negative propagation. 

The phase functions can be regarded as having an overall winding between f = 1 and f = 0, with local phase glitches. 
For given (<r s , ip e i, e s ), the total winding of each perturbed variable across the interval < £ < 1 is a function of harmonic 
number n. The number of winding cycles is not necessarily an integer. Generally the number of cycles increases with n. 
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4-2.3 Illustrative cases 



For the case of (a" s ,i/) e i) = (0.5,0.5) (see Figures B5-B10), the pressure phases ip v and tp c wind positively for the first six 
modes when e s = 0, 1. When e s = 100, <p e may wind negatively in some regions. For n = 1 the winding is negative throughout 
the entire post-shock region. For n > 1 the winding is negative near the lower boundary (f ^ 0.05), but it can be positive or 
negative elsewhere. 

The phase ipQ winds in a monotonically negative sense below £ 0.98 for the fundamental and first overtone. For the 
higher harmonics, the winding of ip^ undergoes one or more abrupt jumps or reversals in narrow ranges of £. There are n — 2 
phase jumps in ip^, and their positions correspond to the distinct density-amplitude nodes. For the harmonics with n > 2, ip^ 
begins with a negative winding near the shock (£ = 1) and remains negatively winding throughout most of the flow, except at 
the jumps. The jumps can be either positive or negative. Descending from the shock, the first jump is positive, and for many 
choices of (<T 3: ijj e i, e s ) and n the second jump is negative. The signs of the further jumps depend on the mode and (a St iJ) e i, e s ). 
For (a s ,tp e i,e s ) = (0.5, 0.5 , 1) t he sequence of jumps is [+,—,—] when n = 5 (e s = 0, 1 in Figure B9) and [+,—,+,—] when 



n — 6 (e s = 0, 1 in Figure B1C). In cases where the cyclotron cooling dominates the positive jump in ip^ occurs very close to 
the shock and the subsequent jumps are all negative and indistinct. 



4-2-4 Phase "discontinuities" 

The sense of a phase jump, e.g. either a modest negative jump in phase or a positive jump by more than half a cycle, 
is essentially distinct because the phase is seen to either increase or decrease asymptotically on either side of the node 
discontinuity, i.e. it is determined by the sign of —dip/d^ in the neighbourhood. In some situations there are critical values 
of e s at which the sense of a phase jump changes from positive to negative. For e s far from the critical values, phase jumps 
are gradual, being spread relatively broadly in £. At the critical e s the phase jump is a discontinuity without a well-defined 
sign. Figure ^ shows an example of the changes in structures of a density phase jump with e s shown near and far from its 
transition value. 



4.2.5 Phases at £ = 

As not all the phases at the lower boundary are directly constrained by boundary conditions, they can only be obtained by 
integrating the hydrodynamic equations. 

The lower-boundary phase of each perturbed variable increases between a mode n and the consecutive mode n + 1, when 
(a s ,tp e i,e a ) is fixed. This inter-mode increment differs slightly between the phases of different hydrodynamic variables. It also 
depends weakly on harmonic number n. (See Figure ^| for examples of lower-boundary phases for one-temperature systems 
with different values of the parameter e a .) 

For a particular mode, there are regular phase relationships between different perturbed variables evaluated at the lower 
boundary. The general trends are as follows: (1) For both one-temperature and two-temperature accretion flows, the density 
and longitudinal velocity profiles are approximately in phase (i.e. ¥?cU=o ~ ]£=())• (2) The phase tp v is approximately 
quarter-phase behind the density and longitudinal velocity phases (ip n \^- Q w ip^\^ =0 — 0.5 « tp T \^ =0 — 0.5). The phasings 
fa — <fib for two hydrodynamic variables a and b, are approximately constant for n > 2. For n = 1, and 2 the phasings conform 
less to the general trend. 



5 BOUNDARY CONDITIONS 

So far we considered a special case of "stationary wall" in which both r and A T are zero at £ = 0. We now consider modifications 
in which the zero-velocity boundary condition for the stationary solution (to = at £ = 0) is retained but the boundary 
value of A T is not necessarily zero. For example, the condition A,r = implies constant (unoscillating) pressure at the lower 
boundary, whilst allowing oscillations of the velocity at the base of the settling flow. Other interesting conditions studied 
include A^ + A T =0 (opposite density and velocity oscillations) and Af = (constant density at lower boundary), (see 
Figure ^). However it is found that the change of boundary conditions has little effect on the frequency sequence <5 : . The mode 
stability S R only changes appreciably in the regime of large e s and high harmonic number n. Some choices with reasonable 
physical interpretations In general most of the choices that we consider yield eigenvalues that are very close to those of the 
conventional case (A r = 0) at £ = 0. Detailed discussions of the effects of boundary conditions on the eigenf unctions can be 
found in Saxton (2001). 

Alternative lower boundary conditions on the stationary variables may be considered in conjunction with alternative 
conditions on the perturbed variables. However these choices would describe systems that are physically very different from 
the white-dwarf accretion problem illustrated in this paper. An exploration of the properties of such systems is beyond the 
scope of this paper, and will be investigated in the future. 
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Figure 3. Phase of perturbed variables evaluated at the lower boundary (£ = 0) for the n = 1, . . . , 8 modes of one-temperature systems. 
Phases are defined in relation to the oscillation of shock position. The left, centre and right columns show cases of the efficiency parameter 
e B = 0.1, 1, 10 respectively. The polar plots show density, longitudinal velocity and pressure phases as diamonds, crosses ans squares. The 
radius is the oscillatory part of the eigenvalue, 8 1 . The bottom row shows equivalent linear plots of phase vs 5 : for the same modes, with 
the same symbols as the polar plots. The phase variables are defined in terms of n, so that a phase difference Aip = 2 means a full cycle. 
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Figure 4. Density eigenfunction of the n = 3 eigenmode of a system with (c s ,^>ei) = (0-2, 0.1), with varying e B . Far from the critical 
value for these parameters (e s ss 6.527), the phase jumps are broad: first row (e s = 10); fourth row (e s = 3). Near the transition, the 
phase jump is abrupt and the amplitude approaches zero at the node: second row (e s = 6.55); third row (e s = 6.5). 



6 LUMINOSITY RESPONSE 

The total power (normalised to p a «g) is given by integrating the cooling function over the whole post-shock structure in the 
stationary solution: 

Jo 7-1 ? J 7- 1 32( 7 -l) V ; 

For the adiabatic index 7 = 5/3, the total power radiated via all processes is L = 1/2, consistent with energy-conservation. 
The contribution of bremsstrahlung cooling is 

L br , = r W-^-ro] 1 dTQ (35) 

Jo 7-1 1 + A cy /Ab t 

The second cooling process (which is assumed to be cyclotron cooling) contributes the difference between this value and the 

total, Z/ C y,0 = L — Z/br, ti- 
lt can be shown that the local luminosity responses are described by complex A-functions that are analogous to the 
eigenfunctions of the hydrodynamic variables. The bremsstrahlung and cyclotron luminosity response functions, normalised 
to e, are 

3 1 

Abi = -A c + - Ac (36) 
and 

A cy = (0.15 -2.5) A C + 2.5A C (37) 

for the conditions of accreting magnetic white dwarfs. For a particular mode the Ab r and A cy profiles describe what could be 
regarded as eigenfunctions of the effect of the oscillations on the cooling emission. 

The amplitude and phase profiles of the modes reveal several regularities. The phases tp^ v and tp cy are both zero at the 
shock, £ = 1 (see Figures meaning that the emission due to both processes locally near the shock oscillates in phase with 
the oscillation of the shock height. The function |Abr| has a minimum in amplitude at or near the shock, and the cyclotron 
function |A cy | has a maximum in the same vicinity. In between the upper and lower boundaries there is no simple relationship 
between the phases of the two functions. Near the lower boundary, both eigenfunctions reach their maximum amplitudes for 
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Figure 5. Eigenplanes for ((Tsji/'ei) = (0.5,0.5) and e s = 0,1,100 with alternative lower boundary conditions. Contours indicate the 
agreement between integrated values of the perturbed variables at given complex <5 and the boundary condition. The top row is for the 
conventional zero-velocity boundary condition A T = 0. The upper-middle row is for X n = 0, meaning constant pressure at the lower 
boundary; the lower-middle row is for A^ + A T = 0, which relates to the strict relation p oc X/v; and the bottom row is for Aj = 0, which 
is a condition of constant density at the boundary. 



any given mode, and these maxima are higher for modes with higher harmonic number n. At the lower boundary, £ = 0, the 
oscillations of the two cooling processes are in antiphase. However, in almost all cases the antinodes of Abr occur at the nodes 
of A cy and the nodes of Abr occur at the antinodes of A cy . There is no obvious relation between the luminosity responses Ab r , 
A cy and the eigenvalues 5. 

Multiplying Ab r and A cy by the respective cooling functions and integra ting over the entire post-shock region yields the 
luminosity responses for a small shock- height perturbation e (see subsection [2.4 ) : 

f 1/4 A br d£ f 1/4 [ 7 (1 - to) - to] A br 
L>br,i = £ / 7^ br 7j — dro = £ / ~ = — Abr<#n> , (38) 



7 - 1 dr ' J 7-I A br + A c 



and 



W = £ ^X cy ^dro = £ Ml - td) - *»»] _ ^ Acydro . (39) 

J 7-1 y dr J Q 7-1 A br + A cy 

L br> i is different for different modes, and so is I/ cy ,i. Moreover, the modes which have high- amplitude Z/ cy ,i oscillations do 
not necessarily have strong oscillations in Lbr,i. 
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Figure 6. Profiles of luminosity perturbed variables, Ab r and A cy , in modes n = 1,2, 3,4 from top to bottom. This choice of system 
parameters, (cr B ,ip e [,e B ) = (0.2,0.1,100), gives strong two-temperature effects and cyclotron cooling dominates the cooling. The left 
column shows the amplitudes for bremsstrahlung luminosity |Ab r | (solid line) and cyclotron |A cy | (dotted line). Central and right columns 
are the phase profiles, ipy >1 and y?cy, with the dotted straight line being a zero- phase reference. 



In Table [j] we show the integrated luminosity responses, scaled according to the amplitude of shock-height oscillation e. 
These quantities are proportional to the relative variations in bremsstrahlung (cyclotron) emission. 

Lbr.i and L cy ,l do not show obvious dependence on <5 R , however they seem to be dependent on the system parameters 
(a s , ipei, e s ). Whether or not |I/br,i|/ibr,o > |icy,i|/^cy,o depends strongly on e s , but is only weakly dependent on cr s and ip e i. 

The complex phases of the values Lbr,i and £ cy ,i, are $br(^) and <3? cy (n) respectively for each mode n, are listed in 
Table ^. Depending on the difference between the phases, the waxing and waning of the emission due to one cooling process 
may follow or lead the other process, or else they may be in phase or antiphase. 

For small e a and for a given mode n, the phases $br(w) and $ cy (n) are nearly constant throughout the (a s ,i/j c i) param- 
eter space. The two-temperature parameters (a s ,tp c i) are almost ineffectual in the small-e s regime. For constant <r s and e s , 
decreasing ip a i causes the phase difference < J > cy (w) — $br(?i) to decrease. This means that if cyclotron luminosity lags then its 
lag increases; else if cyclotron luminosity leads bremsstrahlung luminosity then cyclotron's lead decreases. 

In general, <I?br(l) ~ 0.357T to 0.5-7T and $ cy (l) ~ 0.97T for the fundamental mode in a wide range of system parameters, 
i.e. cyclotron emission oscillation almost always lags bremsstrahlung emission by ~ 0.6n. In the extreme cases in which two- 
temperature effects are so strong that the fundamental mode becomes unstable, e.g. when (<r s , i/i c i, e s ) = (1.0,0.1,100), this 
relation breaks down. The p hase p roperties are more complicated for the overtones because of more complicated winding and 
nodes in the A- functions (see 4.2.2). Therefore no obvious relationships are found between the stability of a mode, the phases 

$br, 



7 TRANSVERSE PERTURBATION 
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7.1 eigen- function profiles 

In the absence of transverse perturbations (re = 0), the transverse velocity eigenfunction X y is zero everywhere. When re > 0, 
the profiles of the other eigenfunctions are modified (see Figures ^, ^). The amplitude |A y | generally has its maximum value 
near the lower boundary, and from that value it declines steeply in f to the value fixed by the boundary condition at the shock. 
Increasing k causes the amplitude \\ y \ to increase in the region near the lower boundary, and the slope d\\ y \/d£ steepens 
throughout the profile. Greater e s makes the slope steeper, when re is fixed. For some « there are one or more local minima 
in the |A y | profile, with the number of minima depending on the harmonic number n. 



The transverse velocity phase ip y generally winds in a negative sense from the shock to the lower boundary (see 4.2.2). 
The total number of turns of cp y increases with re until a threshold is reached. Beyond the threshold there is no winding in ip y 
(see Figures ^-[j, third panel, right column). The lower boundary values of |Af|, |A T |, A„- and |A e | all increase with re beyond 
the re-threshold. 

The |Af| and |A T j profiles have the most distinctive features. However the A<; and A T features disappear when « is 
sufficiently large. These re-dependent properties depend on the harmonic number n and the system parameters (e.g. re ^ 4 for 
(cr s , V>ei, e s ) = (0.5, 0.5, 100) and n = 1, 2, or re <; 8 for (a s , tp ei , e s ) = (0.5, 0.5, 0) and n = 1, 2). 

The amplitude and phase profiles of the total pressure A^ and electron pressure A e are less dependent upon re than the 
density and velocity eigenfunctions are. In the case of (a s ,tp c i) = (0.5,0.5) the eigenfunctions for total pressure and electron 
pressure are not greatly affected by the introduction of a transverse perturbation when re is small. The total pressure and 
electron pressure eigenfunctions are much alike because the amplitude and phase profiles match at both the upper and lower 
boundaries. Increasing re causes decrease of the lower boundary pressure phase, y B |f=o = Vvrl^o- The phases <p-n and ip e have 
similar profiles for small-e s flows, but the electron pressure develops phase jumps and associated node-like amplitude features 
when e s is large. (For a more detailed discussion of the transverse perturbation, see Saxton 1999.) 
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Table 1. Amplitudes of the oscillations of the total emission in bremsstrahlung and cyclotron radiation, relative to these processes total 
luminosities in the stationary solution. For each set of the system parameters, the left column indicates whether the mode is stable (— ) 
or unstable (+) in the small-amplitude analysis (see complete eigenvalue results in Saxton & Wu 1999). The middle and right columns 
are the e-normalised relative amplitudes I/br,l /eibr.o an d |-^cy,i| /e£cy,o respectively. 
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7.2 Luminosity response 

Figures |lo| - |r"i| show the luminosity responses as functions of k for n — 1 and 2. The integrated luminosity amplitudes |£t>r,i| 
and | icy, 1 1 have minima in « where there are abrupt changes of the phases <E>br and $ C y respectively. The number of K-minima 
is determined by n. For n = 1 and (er s , ip e i) ~ (0.5, 0.5) the |-E<br,i| minima occur at about k ~ 1.6, 1.0, 0.6 for e s = 0, 1, 100 
and the |£ C y,i| minima coincide at approximately the same K-values. For n — 2 there are at most two minima of |Lb r ,i|- For 
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Table 2. Total brcmsstrahlung and cyclotron emission phases, and $ cy , expressed as multiples of it, for the modes n = 1 . . . 6 from 
top to bottom, under given (<r s , ip c i, £ s) conditions. In this convention the phase of the shock-position oscillation is defined as <I> S = 0. In 
each set, the third column described the difference between the phases of integrated cyclotron- and brcmsstrahlung-luminosity oscillations, 
*cy(i) — *br( n )- Negative values indicate the cyclotron luminosity oscillation following the brcmsstrahlung luminosity oscillation; and 
positive values indicate cyclotron luminosity leading bremsstrahlung luminosity. 
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small e s the minimum corresponding to a larger K-value becomes an inflection, e.g. for e s = the actual amplitude minimum 
of |ibr,i | is at k ~ 1.8 and the inflection point is at k « 3.4. 

Both amplitudes, |£br,i| and |L cy ,i| are slowly varying in k when the value of n is below k», where n*(n, e s ) ~ ^(2n — 
l)/t e (e s ) and K t depends on the system parameters. For (a s ,tp e i) = (0.5,0.5), n e ~ 1.0 for e s = 0. It reduces gradually as 
the cooling efficiency increases, and k € ~ 0.3 at e s = 100. For wavenumbers n > k», both |Lbr,i| and |L cy ,i| increase rapidly 
with k. This is a consequence of the increasing amplitudes of the hydrodynamic variables' A-functions as k increases. Because 
bremsstrahlung cooling is most efficient in regions near the lower boundary, |Lbr,i| tends to rise more steeply with n than 
|L cy ,i| does. 

The integrated luminosity phases $br and <3> cy wind with k. On top of these winding trends, there are phase jumps where 
the respective £br,i|, £c y ,ij amplitudes reach minima. Away from the minima, both of the phases generally decrease when k 
increases. Generally the winding of both phases in k is more rapid when e s is large, however for sufficiently large e s there are 
modes where 3>br overtakes <E> cy in its variation with n, e.g. n — 2 with (<r s , ip e i, e s ) = (0.5, 0.5, 100). The prevailing winding of 
the bremsstrahlung phase <3> br is usually more sensitive to k than the <3> cy is, i.e. — d$ hT /dK tends to be greater than — d$ cy /dft. 

In summary, the presence of transverse perturbations may significantly alter the instability of a mode and modify the 
luminosity response. 
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Figure 8. Eigcnfunction profiles in the presence of transverse perturbation, for density, longitudinal velocity, transverse velocity, total 
pressure and electron pressure (from top to bottom), for the n = 1 mode with (<Ts,i/ , ei I ^s) = (0.5,0.5,1). Profiles for wavenumbers 
K = 0, 1, 4 are plotted in dotted, dashed and solid curves respectively. The left column shows logarithms of the amplitudes, and the right 
column shows the phases. Phases are multiples of w. The dotted line in the phase plots gives reference for zero phase, which is defined 
to be the phase of the oscillation of shock height (marked in grey). 



8 CONCLUSIONS 

We have presented a general formulation for the linear analysis of two-temperature radiative shocks with multiple cooling 
processes. The formulation recovers the restrictive cases in the previous studies such as the one-temperature flows with a 
single cooling function (Chevalier & Imamura 1982), the one-temperature flows with multiple cooling processes (Saxton et al. 
1998) and the two-temperature flows with a single cooling function (Imamura et al 1996). We have applied the formulation to 
mCVs and investigated the hydrodynamic and emission properties of the time-dependent post-shock accretion flows in these 
systems. Our finding are summarised as follows. 

The amplitude profiles of A-eigenfunctions show local minima and maxima, which we identify as nodes and antinodes. 
The nodes and antinodes are prominent only in the eigenfunction profiles of density and electron-pressure. The eigenfunctions 
for the fundamental and first overtone are more similar to each other than any of the higher overtones. The eigenfunctions 
for higher-order modes have more nodes. 

The phase profiles of the eigenfunctions describing particular hydrodynamic variables circulate about the complex plane 
as £ varies from the shock down to the lower boundary. This circulation can be positive or negative overall, or there may be 
reversals of the winding sense between distinct zones. The abrupt jumps in the phase profiles always coincide the nodes in 
the A-eigenfunction profiles. 

The luminosity responses of cyclotron and bremsstrahlung are determined by the A-eigenfunctions. There is no obvious 
general relationship between the amplitude or phase luminosity responses and the stability properties of the flow. There are 
situations in which a mode is unstable but the amplitude of oscillations are larger in the cyclotron luminosity than in the 
bremsstrahlung luminosity. 

For the same stationary condition (r = at £ = 0), all out choices of perturbed boundary conditions do not show significant 
differences in the stability properties. We therefore conclude that the stability properties of the flow with stationary-wall 
boundary is mainly determined by the energy-transport processes. 
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Figure 9. Same as Figure |i] but for n = 2. 



The presence of a transverse perturbation modifies the eigenfunction profiles of all the hydrodynamic variables. The 
profiles of electron-pressure and total-pressure eigenfunctions are less affected in comparison with the other eigenfunction 
profiles. In some range of transverse wavenumber the density and longitudinal- velocity eigenfunction profiles develop extra 
node features. When n is large enough (Ss 3 for n = 1,2), the amplitudes of all eigenfunctions become large near the 
lower boundary. The amplitudes, however, decrease as the height £ increases. For some values of the transverse wavenumber 
(1 ?S K ^ 3), a mode which is stable in the absence of the transverse perturbation can become unstable. However, when k is 
very large, the mode is stablised. The phase difference between the oscillations in the bremsstrahlung and cyclotron luminosity 
are also modified in the presence of transverse perturbations. 
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Figure 10. Integrated luminosity response amplitudes and phases for n = 1 mode in the presence of transverse perturbations, for the 
parameters (cr s ,t/> c i) = (0.5,0.5) and e s = 0,1,100 from left to right. The curves corresponding to bremsstrahlung cooling are marked 
with +; the curves corresponding to cyclotron cooling are marked with X. 
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Figure 11. Same as Figure flu but for the n = 2 mode. 



APPENDIX A: REDUCTION TO RESTRICTED SYSTEMS 

The matrix equation describing the perturbation is 
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For systems with a purely longitudinal perturbation, k — and the third row and third column of the matrix can be 
eliminated. The \ y terms vanish from the remaining F functions, and the equation for \ y decouples from the rest of the 
perturbed variables is simply (lnA H )' = S/to — (lnro)' or otherwise X' y = X y (S — Tq)/to- The two-temperature system with 
purely longitudinal perturbations is described by this reduced matrix equation, with corresponding F functions that omit all 
terms of X y : 
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In the single-temperature limit, i/j e i becomes large and the electron and ion pressures both equal half of the total pressure, 
i.e. we have 2ir E — > 7ro = 1 — To, and A c — > A^ throughout the entire post-shock flow. Moreover, F$ — > hF&. Then we can 
eliminate the fifth row of (|Al|), yielding 



Stability analyses of two-temperature radiative shocks 21 



d 
dr 







1 -1 





Xr 


1 


-j-kq/to 1 





\y 


~~ A 





-(77T - 


Att 




7 -7 






1/ro 




" Fx 


-1/t 




F 2 







F 3 






Fi 



(A3) 



Reducing the system to a one-temperature form with purely longitudinal perturbations leaves only three non-trivial perturbed 
variables and hence a 3 x 3 coefficient matrix 
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which is equivalent to that in Saxton et al. (1997) and Saxton et al. (1998). 
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APPENDIX B: EIGENFUNCTION PROFILES 





Figure B2. Same as Figure Bl but with n = 2 and (cr s ,ip e i) = (0.2,0.1). 
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Figure B4. Same as Figure 



Bl 



but with n = 4 and (<r s ,i/> ci ) = (0.2,0.1). 
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Figure B6. Same as Figure 



Bl 



but with n = 2 and ((r s , ipel) = (0.5, 0.5). 
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Figure B8. Same as Figure 



Bl 



but with n = 4 and ((r s , ipel) = (0.5, 0.5). 
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Figure BIO. 



Same as Figure FS]I but with n = 



6 and (cr s , ip ci ) 



= (0.5,0.5). 



